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ABSTRACT 



Three approximation methods for generating outcomes on 
Poisson random variables are discussed. A comparison is made 
to determine which method requires the least computer execu- 
tion time and to determine which is the most robust approxima- 
tion. Results of the comparison study suggest the method to 
choose for the generating procedure depends on the mean value 
of Poisson random variable which is being generated. 
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I. INTRODUCTION 



It is frequently desired to generate Poisson random 
variates in simulations. There are standard exact methods 
for doing this; the problem arises when a computer is used 
to generate the Poisson random number which has a large mean. 

For example, generating one random number such as 105 from a 
Poisson distribution with mean 100 needs at least 105 calls 
to a pseudo-random number (uniform (0,1)) generator. Compu- 
ter time requirements become important cost factors when 
considering various methods for generating random numbers. 

The objective of this paper is to examine several ap- 
proximated ways of generating Poisson random variables and 
to determine the method which gives minimum execution time 
and small mean absolute deviation according to the Poisson 
mean value. The mean value is the only parameter in the 
distribution. Comparison statistics to determine the best 
approximation to the Poisson distribution are the cumulative 
probability, mean absolute deviation and Kolmogorov-Smirnov 
test. Finally a composite generating procedure according to 
the mean value is suggested. 

The following notation is used in this study, 
denotes a uniform (0,1) random variable; 

N denotes a Poisson random variable where mean is m; 

Z denotes a random variable from a unit normal distribution. 
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II. GENERATION OF POISSON DISTRIBUTED VARIATES 



A. THE POISSON DISTRIBUTION 

A random variable N with integer values has a Poisson 
distribution if 

-m n 

prob ' {N = m} = e ^ n = 0,1,2... 

1 1 • 

In order to generate a Poisson random number N from a Poisson 
distribution with mean m, the following algorithm is pre- 
sented. It is. the standard exact method for generating these 
deviates . 

Let U^, i = 1,2,..., be independent uniform (0,1) random 
variates. The Poisson variate, N, is computed as: 



N 




if U. < e m 

i - 

k k+1 

if n U. > e > H U. 

. , i - . , i 

i=l i=l 



where N is non distributed as a Poisson with mean m, i.e., 
Prob (N=n) = e m m n /n! . 

Equivalently, since the logarithm is a monotone transforma- 
tion, we have 



N 



0 if ln(U^) < ln(e m ) 



wn 



k 



k k 

if ln( n U ) = E ln(U jL ) 
i=l i=l 

k+1 k+1 

-m > E ln(U. ) = In ( II 
1=1 1 1=1 



> ln(e~ m ) 
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Changing the signs and direction of the inequalities we 
have : 



N 



or letting = 



N 



0 if -ln(U i ) > m 

k k+1 

Ik if E ln(U. ) < m < E 

\ i=l 1 i-1 

- ln(U ± ) 



r 

; 0 if E. > m 

l - 




if 



k 

I E. 
1 



k+1 

< m < E E. 
- i-1 1 



ln(U ± ) 



where the E^ are exponentially distributed variates with 
mean 1. 

If n multiplications of uniform (0,1) random numbers is 
strictly greater than e -m and if n+1 multiplication of uniform 
(0,1) random number is equal and less than e -m then n is the 
Poisson random number. Generally generating one random num- 
ber from Poisson process with parameter m requires on the 
average m+1 uniform (0,1) random numbers. This is because 
the number generated is n+1. When m is large it is clear 
that generating Poisson random numbers with the above method, 
although it is exact takes a lot of computer time and this 
method may be uneconomical. In addition, the large number 
of multiplications can produce serious precision problems 
on a digital computer. 
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B. APPROXIMATIONS FOR POISSON VARIATES 



1. Normal Approximation 

In a Poisson process with parameter X it is necessary 

to generate random variables from the Poisson distribution 

with parameter (mean) m. Now look at counts in (0,x), where 

x satisfies Xx = m. The central limit theorem says that as 

m goes to positive infinity, (or when x gives the infinity in 

Poisson process with fixed X), then N, which has mean m and 

u 

variance m, is such that N+^-m/m is approximately distributed 
as a unit normal random variable. Denote a random variable 
from a unit normal distribution by Z . So N is distributed 
approximately as m^Z + m. In order to generate Poisson ran- 
dom numbers from the normal distribution, first generate Z; 
then let 



N 



1 0 if m^Z + m - 0.5 < 1 

| m^Z + m - 0 . 5 I otherwise 



where Jl§J denotes the greatest integer less than or equal to 
a; also known as the "floor" of a. N is then the approximated 
Poisson random variate. 

2 . Square Root -Transformation of Poisson Distribution 

If N is a Poisson random variable with mean m then 

Y = /N + 3/8 is approximately distributed as a normal dis- 

u 

tribution with mean m and variance h. This result is due to 
Bartlett [Ref. 9]« 



11 



This method is derived as follows: let Y = /N + C ~ N(y, 0 2 ) 

where C is a non-negative constant. Let t = N - m and 
m* = m + C. Define coefficients for s = 1,2,3* •••by 

a = (_i) s+1 l*(-l)*(-3)***(-2s+3) 

s K ' 0 s , 

2 • s : 

Then for any t > - m' we have a Taylor series expansion. 



If t > 0, 
is bounded 
of t are y 



Y = /m 7 " {1+A- - A 0 (^i-) 2 +• • * + (-l) s A ,(^r) S+1 } 

1 m f 2 m f s-1 m 1 



+ V 



we see at once | R | < A t s /(m') s ^ converges and 

s s 

(P. J. Anscombe (4)). We note now that the moments 
! = 0, y 2 = m, y 3 = m, yi* = 3m 2 + m,..., which give 



var ( V . k(1 + 3gl* 32C ‘ + 17 >, 

;o that when C = 3/8 Var (Y) ~ \ (1 + l/l6m 2 ). Also 



E(Y) 



Sm+C - 



1_ 24C - 7 

8m* 128m 372 



Let XNR = /N + 3/8. Then XNR is approximately normally dis- 
tributed with mean /m and variance h. 



„ _ XNR - /m + 3/8 

h 

XNR = | + /m + 3/8 , 
/N + 3/8 = | + /m + 3/8 , 
thus set 
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N ■ 



j(f + /m+37FT 2 - 3/8 



if (| + /m+3/b) 2 - 3/8 < 1 



otherwise 



N is then the approximated Poisson variate. 

We now need to calculate the probability distribution 
of N obtained in this way from the square root transformation. 
We want the probability that 



n - 1 + 3/8 < XNR 2 < n + 3/8 



if we divide by the variance % we get 

M(n - 1 + 3/8) < XNR 2 4 < 4(n + 3/8). 

Note that 4*XNR 2 is distributed with a non-central X* distri- 
bution. The non-central Xi density is 



f (x) = 
x' 



-Jg(x+y 2 /a 2 ) - y 



2/x /2i 



[e 






+ e a 



H/x 



]. 



Thus if y=m^ } o 2 =(%) 2 s then 



and 



f x (x) 



e -ia(x+niA) + e (M/H) ] 

2/x STv 



prob(x=N) = 



l 



N+3/8 

(h ) 2 



f (x)dx - 
x ' 




N+3/8-1 

Lh ) 2 



f (x)dx 
x v ' 



This allows us to evaluate directly how well the 
distribution of N approximates the distribution of a Poisson 
variate with parameter m. 
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Note that since in the LLRANDOM package it takes the 
same amount of time to generate 5 uniform random variables 
as it takes to generate a normal random variable, the pro- 
cedures will be competitive timewise once m is much greater 
than 5 . 



3 . Cube Root Transformation of Poisson Distribution 
If N is a Poisson random variable with mean m then 
Y = ftj-1/24 is approximately distributed as a normal distri- 

"I/O 

bution with mean m ' ^ and variance 1/9 3 /m when N > 1. This 
comes from the following derivation which is essentially the 
same procedure as for the square root transformation. Sup- 

3 



pose 
m 1 = 



Y = /N+C is distributed as N(y,a 2 ) 
m + C. Define coefficients for s 

a = (-1 ) S+1 l'(-2)(-5)(-8)'* y4 



. Let t 

= 1,2,3, 

-3s) 

T # 



= N - m and 

. . by 



For any t > - m* we have the Taylor series expansion 



Y - for U + ai s* - a 2 (jjtt) 2 + > 3 ( 5 *) * - 



t \ s— 1 ■ 



• + <-« vi tfy > + 



R 



if t > 0, we see at once that | R t /(m') S converges. 

S s s 

Therefore , 



f(s) (1 + e ST ) t s 
R s = JT (5) 0 1 9 < 1 



It | < 



m ’ 



R (m 1 ) 1/3 = (1 + -r ) 1/3 
s m' ' 



- 11 * a l F -•••+(-!) 



/ t \ s 1 
a s -l ( E> 
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R (m') 1/3 » 1+1 t i 

= Z (-1) 1 a. (— r) 1 converges and is 

t s i=s 1 m 



bounded . 



We note now that the moments of t are 1^=0, y 2 =m, 
y 3 =m, y 4 = 3m +m, — } giving us 



and 



E(Y) = 3 V^n+C" - +... 

m ^ 



Var(Y) = ^/m 1 ") 2 { xm “ 



1 2m 2 +m 
81 (m')“ 



9 3 Sm 



i ( JLc + i ) 



27 “ ■ 152 



If C « 



then 



^2A 
Var(Y) = 



9 3 Sm 



Let 



YNR = 3 /N - 1/24 N > 1 



Z = 



YNR - 3 /m + 1/24 



3^ 



YNR = 3 /m + 1/24 +(l/3 6 >/m) Z 



3 /N-l/24 = Vm + 1/24 + (l/3 6 /fn) Z. 



Thus set 



0 



if ( 3 /m+l/2^l+| 6 /m Z) 3 +^ < 1 



N = 



| ( 3 /m+l/2*f + 3 6 Sm Z ) 3 + 



otherwise 
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N is now the approximated Poisson variate. 
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III. EVALUATION OP THE METHODS 



Generally, generating Poisson random number from the 
exact method is known to take a long execution time, since 
one generated random number "N” requires on the average N 
multiplications of uniform (0,1) random numbers. Therefore, 
, generating Poisson variables with the exact method (which 
gives the best accuracy) is good for a small mean, m, while 
the approximation methods, which take shorter execution time 
but with less accuracy, should be used for large m. Here we 
need a trade-off between execution, time and accuracy to 
choose the generation method according to the mean value of 
the Poisson distribution. 

The comparison statistics show how closely the methods 
approximate the original Poisson distribution. In the 
Kolmogorov-Smirnov test, all approximations are accepted at 
significance level a = .05. Prom the comparison of the em- 
pirical probability distributions and the mean absolute 
deviations from the exact distribution, the optimal genera- 
tion procedure based on the mean value, m, is as follows: 

Method Mean (mj 

1. exact Poisson distribution if 0 < m < 20 

2. square root transformation .if 20 < m < 100 

3. normal approximation if m > 100. 

The cube root transformation should not be adopted be- 
cause it is far less accurate compared to the square root 
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transformation and the normal approximation. The following 
tables and figures show the comparison statistics of the 
three approximations versus the exact Poisson distribution. 

Table I assesses the accuracy of the normal, square root, 
and cube root approximations to the exact Poisson distribu- 
tion at selected points. Figure 1 illustrates graphs of the 
empirical cumulative distribution functions of the three 
approximations to that of an exact Poisson distribution with 
mean 10. Figure 2 shows graphs of the mean absolute devia- 
tions of the three approximations from the exact Poisson 
cumulative distribution function for various mean values. 

The mean absolute deviation is defined as: 



where 




1 

k-1 



k 

Z 

i=l 




P.) 2 } 



Pj^ is the CDF of the approximating distribution; 

P. is the CDF of the exact Poisson distribution; and 
i 

k is the sample size. 

Table II is a summary of a comparison of the empirical 
distributions produced by the three approximation methods to 
the exact Poisson distribution by means of the one sample 
Kolmogorov-Smirnov test. The null hypothesis is that the 
approximated distributions are Poisson against the alterna- 
tive that they are not Poisson. 

Lastly, Table III and Figure 3 analyze the sensitivity of 
the square root transformation to the constant C. In the 
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derivation of the square root transformation, C was chosen 
as 3/8. Different constants were used in order to find the 
most robust constant to use, i.e., the constant which yielded 
the smallest mean absolute deviation from the exact Poisson. 
The value C = 13/18 was found to be the most robust. Note 
that this value of C was used in the approximation when 
making the comparisons with the other methods. 
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Table I. Accuracy of the Normal, Square Root, and Cube Root 
Approximations to the Exact Poisson Distribution. 









p < z (m) i "> 


* 


Poisson 
mean, m 


Observed 
Value, n 




Normal 


Square Root 
Transformation 


Cube Root 


Exact 


Approx 


0=13/18 


Transformation 




2' 


0.08422 


0.07301 


0.07719 


0.10491 


5 


3 


0.14037 


0.11939 


0.14051 


0.18491 




6 


0 . 14622 


0.16036 


0.14614 


0.16681 




10 


0.04861 


0.04485 


0.04704 


0.09045 


15 


14 


0.10244 


0.09937 


0.10288 


0.05512 




18 


0.07062 


0.07622 


0.07073 


0.06348 




l4 


0.03874 


0.03633 


0.03746 


0.08819 


20 


19 


0.08884 


0.08683 


0.08891 


0.08726 




22 


0.07692 


0.08050 


0.07672 


0.05755 




36 


0.05394 


0.05161 


0.05409 


0.02552 


40 


39 


0.06296 


0.06223 


0.06307 


0.06336 




42 


0.05850 


0.05995 


0.05841 


0.05675 




55 


0 .03960 


0.03910 


0.03963 


0.00662 1 


61 I 


59 


0.05019 


0.04999 


0.05020 


0.04973 




82 


0.04407 


0 .04402 


0.04402 


0.03289 


82 


87 


0.03686 


0.03649 


0.03683 


0.03801 1 




89 


0.03165 


0.03243 


0.03149 


0.03991 




94 


0.03775 


0.03805 


0.03759 


0.0255 


90 


97 


0.03H2 


0.03182 


0.03098 


0.0435 




93 


0.03223 


0.03241 


0.03217 


0.0410 


100 


99 


0.03997 


0.03980 


0.03994 


0.03605 




112 


0.02881 


0.02847 


0.02866 


0.03627 


120 


116 


0.03477 


0.03438 


0.03462 


0.00130 
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Observed n | 

Figure 1. Comparison of Empirical Cumulative Distribution Functions for Poisson Distribution 
with Mean 10. 



Z O’ 
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Figure 2. Mean Absolute Deviations of the Three Approximations. 



Root Cube Root | Critical 
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Table II. 

One Sample Kolmogorov-Smirnov Test. 



c 


/l k 

Mean Absolute Deviation!/ — E (P! - P.) 2 

fk- 1 1= i i i 


Poisson Mean = 20 


Poisson Mean = 80 


3/8 


0.019 


0.009 


5/9 


0.016 


0.006 


7/12 


0.014 


0.005 


11/18 


0.011 


0.005 


23/36 


0.009 


0.004 


8/12 


0.007 


. 0.003 


*13/18 


0.006 


0.002 


9/12 


0.007 


0.003 


5/6 


0.011 


0.004 



* optimal constant for C. 



Table III. 

Sensitivity of 'the Square Root Transformation to the Constant 
C for Poisson Distributions with Means m = 20 and m = 80. 
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mean absolute deviation 




Figure 3. Plot of Sensitivity Test Given in Table III. 
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IV. THE METHODS OF DIETER AND AHRENS 



During the writing of this thesis, a pre-publication 
chapter from a book by U. Dieter and J. H. Ahrens was re- 
ceived. This chapter deals with the generation of non-uni- 
form variates and contains a section on the generation of 
Poisson variates. Some of the methods which’ they described 
will be outlined here. They have not been programmed or 
tested against the composite method of this thesis. 

A. UNIT MEAN METHODS 

It is well-known that the sum of independent Poisson 
deviates with mean is Poisson distributed with mean 

y i+y 2 +. . • Hence the following algorithm may be considered. 

1. Split the mean m of the Poisson distribution into 
m’-«-[m] and f-*-m-m’<l. (Hence [m] means take the integer part 
of m and f is then the fractional part.) 

2. Take m* samples from the Poisson distribution of 
mean 1 using some suitable method. Let s be the sum of these 
samples . 

3. Take a further sample n from the Poisson (1) distri- 
bution. If n = 0 deliver k«-s. Otherwise, generate Ui,U 2 ...U 
and count the number, j, of U. for which U. < f, and deliver 
k-<-s+j . 

The validity of step 3 is proved as follows : the probability 
of observing j samples less than or equal to f (once n - j 
samples greater than f) is 



26 



( j ) r j (i-f) n " J . 

Since n follows a Poisson distribution with mean 1 this 
yields 



p j - 



CO 

z 

n-j 



-1 



ni 



( . ) f J (1-f) 



n-j 



-1 



= Z 



n-J 



f^(l-f) n ^ 



e " (l-f) n _ e 1 f J ’ 1-f _ e f f J ’ 

j I n _ j =Q (n-j)! j : 6 3 l ' 

The unit-mean method becomes a suitable algorithm for 
the case when the exact method should not be used for gen- 
erating Poisson random variates a that is 3 when the mean is 
large (m > 200). This method will then become slightly faster 
than the exact method, but still needs to call too many uni- 
form deviates since it is basically the same as the exact method 
except for splitting the mean, compared to the square root 
transformation approximation method, the unit mean method may 
be uneconomical. 

B. THE CENTER TRIANGLE METHOD 

This method is suitable in the case of variable means m 
of the Poisson distribution. The mean is split into m* = [m] , 
the integer part of m, and the remaining fraction (Hf^-m'cl. 

The algorithm uses a mixture of two distributions. The first 
distribution has an easy sampling method and is employed most 
of the time. It is based on the following isosceles triangle. 
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g(x) = / - - Ix-m' | if m'-s < x < m'+s 

(B. 1) J 3 

*0 elsewhere. 

The constant s is chosen as 
(B. 2) s = 2.2160358671665 - 0.78125. 



The area under g(x) is b . If b is small enough then the 
inequalities 



(B. 3) 



t . 

1 



■f 



g(x)dx < P. (i = 0,1,2,...) 



- 1 



may hold for all i. The maximum b for which this is true was 

found numerically for all ^ < m < 200. The choice of the 

constant a = 2.2160358671665 in (B.2) is well motivated. The 

vertices of the maximum isosceles triangle that can be placed 

wholly inside the standard normal curve are situated at (±a,0) 

and (0,1//2 tF). As m' goes to infinity, b approaches the 

max 

area of this triangle which is .884070^0229876. The correc- 
tion -0.78125 = -25/32 in (B.2) improves the value of b 

max 

for small m. [7, page 11-15, table.] Here the variable mean 
m is split into m = m’tf as in the unit mean method. If 
m < 8 the exact method is applied. If m > 9 , the triangular 
method is used with the triangular probabilities t^ = 27/32. 
This center triangle method needs a long program including 
several tables which requires a large memory space and suf- 
fers from the roundoff errors ' of b . 
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This method proves a very efficient algorithm based on 
a simple idea. However, it is somewhat messy to describe 
the program and has little mathematical appeal compared to 
the square root transformation method. This method also 
needs a long execution time as does the exact method. For 
example, to compute one Poisson variate from a distribution 
with mean m = 100, the program would require 1886 ys [7» 
page 11-19] • 

C. THE GAMMA METHOD 

In order to obtain a sample k from the Poisson distribu- 
tion with variable mean m, select a positive integer or 
(typically n is a little smaller than m) . Then, take a sam- 
ple x from a Gamma distribution with parameter n. 

Case (1) if x > m, return a sample k from the binomial 
distribution with parameters n-1, m/x. 

Case (2) if x < m, take a sample j from the Poisson dis- 
tribution of mean m-x and return k-«-n+J . 

The sample x simulates the n-th event (arrival) in a 
Poisson process of rate 1. If x > m, then there are n-1 
arrivals in the interval (0,x), and each of these has a 
probability of m/x of being below m (Case (1)). If x < m, 
then the n simulated arrivals are all before m and the sample 
j indicates the additional events between x and m (Case (2)). 

A formal proof of the procedure runs as follows. In the 
first case one has m < x < 03 and 
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p 



k 




, , , , -x n-1 

n-l N ,m s k /n m^n-l-k e x 

k ; V U "x ; - r(n) 



is the probability of obtaining k from the binomial distribu- 
tion (n-1, m/x) summed over the Gamma (n)-distributed values 
of x above m. The expression transforms into 



P, = 



(n-1) i 
k ! (n-l-k) ! 



m K (x-m) 

n-1 



n-l-k 



r' x x "- 1 

(n-1) ! 



dx 



m 



k! (n-l-k) ! 



(x-m) n 1 k e X dx 



m 



k -m 
m e 

k! (n-l-k) ! 



t"- 1 ^ e-* dt 



where t = x - m. Notice that f(n-k) is (n-l-k)!. Hence 



k -m 

P = m e 

k k! 



as required. In the second case one has 0 < x < m, and 



P 



k 




(m-x) 



k-n -(m-x) 

e 

(k-n) ! 




x 



n-1 



Ut J 



dx 



o. 

is the probability of obtaining j = k-n from the Poisson 
distribution of mean m-x summed over. the Gamma (undistri- 
buted values of x between 0 and m. This expression trans- 
forms into 
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e~ m f 

-n) ! (n-1 ) T J 



x n ■*" (m-x)^ n dx. 



Introducing t = x/m yields 



-m n-1 k-n 
p _ e m m m 

k (n-1) I (k-n) ! 



S 



t "- 1 d-t) k - n at. 



The integral is a Beta function 3(n,k-n+l) 



p _ e m m k r(n) r(k-n+l) e m m k 

k (n-1)! (k-n)! f(k+l) k ! 



as before. The following algorithm is considered. 

1. Initialize k*-0 3 w-*-m 

2. If w > C (C = 24 was used as the mean cut-off point) 
go to 6 

3. (Start Case (2)). Set p«-l and calculate b«-e -w 

4. Generate U and set p«-pU. If p<b deliver k 

5. Increase k-*-k+l and go to 4 

6. Set n«-[dw] where d = 7/8. Take a sample x from the 
Gamma-(n) distribution. If x>w go to 8. 



7. 


Set k k+w. 


w w-x and go to 2 


8. 


(Start Case 


( 1) ) . Set p«-w/x 


9. 


Generate U. 


If U<p increase k«-k+l 


10. 


Set n«-n-l. 


If n>l go to 9 


• 

i — 1 
i — 1 


Deliver k. 





The performance of the algorithm depends on the cut-off 
point C in Step 2. Step 3> Step 4 and Step 5 are exactly the 
same as the steps of exact methods. This gamma method has a 
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complex sub algorithm for the binomial distribution In Steps 
8-10 and requires memory space for the gamma distribution. 
This method may be efficient in cases with extremely large 
means, .say m = 1000, but this method requires a rather long 
execution time as does the exact method. For example, to 
generate a sample from a Poisson distribution with mean 100 
takes 1175 ys, see [7, page 11-23]- 



/ 
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